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ABSTRACT 

We present a new code for radiation transport around Kerr black holes, includ- 
ing arbitrary emission and absorption mechanisms, as well as electron scattering 
and polarization. The code is particularly useful for analyzing accretion flows 
made up of optically thick disks and optically thin coronae. We give a detailed 
description of the methods employed in the code, and also present results from 
a number of numerical tests to assess its accuracy and convergence. 

Subject headings: black hole physics - accretion disks - X-rays:binaries 



1. INTRODUCTION 

Since the beginning of X-ray astronomy over 50 years ago, there has been steadily 
growing interest in relativistic radiation transport. Because of the high energies of both 
photons and electrons relevant to these astrophysical sources, special relativistic effects must 
be included in most particle interactions. And because the central engines of so many X-ray 
sources are compact objects such as neutron stars and black holes, full general relativity 
must also be incorporated into any physically realistic code. 
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Here we present in detail for the first time the fully relativistic Monte Carlo radiation 
transport code Pandurata0. While this is the first formal description of the code in the liter- 
ature, it has been under development for r nany years and has already been used in numerou s 
publications dSchnittman k Krohkl 120091 . boid iNoble et al.lboili ISchnittman et aljl2012[ ). 
Pandurata shares features with numerous existing codes in the literature, but we believe its 
combination of full general relativity, wide range of emission and absorption processes, scat- 
tering, polarization, and optically thin/thick capabilities make it uniquely valuable in the 
rapidly evolving field of black hole astrophysics. Most recently, in ISchnittman et al.l (120121 ) 
we have demonstrated how the code may be used to take a major step towards bridging the 
gap between global magneto-hydrodynamics (MHD) simulations and real X-ray observations 
of accreting black holes. 

The literature of radiation transport in astrophysics is extremely broad and includes 
scores of different techniques and applications. It would be a futile endeavor to attempt 
to give a comprehensive summary of work here. Rather, we will simply highlight a few re- 
cent contributions that are most relevant to the applications of interest, namely photon 
transport around Kerr black holes. By far the most common approach has been ray- 
tracing geodesic paths backwards from a distant observer to the accretion fiow, calcu- 
lating a transfer function of some sort, and coupling this to some model for emission 
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A smaller number of codes have been written with the emitter-to-observer approach, 
which may be more physically intuitive, but is a lmost alway s more computa tionally in- 
tensive. With the exception of iLaor et al.l ( 1l990l ): iLaorl (Il99ll ): iKojimal ( 1l99ll ). which use 
uniform sampli ng of emission angles , these codes are generally Monte Carlo in nature, such 
as grmonty by iDolence et al.l ( 120091 ). and the present work. As we will show below, par- 
ticularly when electron scattering is included, the emitter-to-observer paradigm is almost 
essential for capturing the most relevant physics of the problem. 

Another feature that is relatively uncommon in these ray-tracing codes, but of increasin g 
interest in the high energy community, is polarization. It is included in lAgol fc KrolikI ( 120001 ): 



^The name Psindurata comes from Coelogyne pandurata, a species of black orchid native to Borneo. Much 
of the core radiatio n transport is derived from t he co de Buttercup, developed by Schnittman for inertial 
fusion applications ( Schnittman fc CraxtonlE996 . 2000 ) at the University of Rochester's Laboratory for Laser 
Energetics (LLE). In holding with LLE's long tradition of naming codes after flowers, the name Pandurata 
was chosen to represent the joint heritage of black holes and laboratory radiation hydrodynamics 
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Dovciak et all (120081 ) .iHuang et al.l (120091 ) : IShcherbakov fc Huand f l201lh : lHuang fc Shcherbakov 
( 120111 ) ; iMarin et al.l ( 120121 ) , although o ften only for yacuu n i transport and no t including scat- 
tering. Disk polarization is treated by iLaor et al.l ( 119901 ): iMatt et al.l (119931 ): iDovciak et al. 



(l2008l) for bo t h Sch warzschild and Kerr black holes, but neglecting electr on scattering 



Dovciak et al.l (120111 ) includes illumination from a source above the disk, while iDovciak et al. 



( 120081 ) includes a cold plane-parallel atmosphere above the disk, geometrically thin with vary- 
ing optical depth. A small number of ray-traci ng codes also allow for non-standard black 
hole metrics as a way of testing general relativity. iKrawczynskil ( 12012| ) fo llows the emitter-to- 
observer paradigm for calculating polarized flux from a thermal disk, and lPsaltis fc Johannsen 
( 120121 ) describes an observer-to-emitter framework that can be applied to a large nu m ber of 



space -time tests such as timing, spectra, and imaging ( iJohannsen &: PsaltisI l2010al Jbl. 12011 
20121 ). 



The body of literature including detailed scattering and pola rization is genera. 



stricted to fiat spacetime, and often only the most simple geometries (|Connors fc Stark 



Connors et al. 



1980 



Poutanen fc Svensson 



Sunvaev fc Titarchuklll985l : iHaardt fc Maraschilll993l : iHaardt et al. 



y re- 



977 



1994; 



19961 ) . Here we attempt to synthesize the strengths of all these various 



codes into a single flexible radiation transport tool for analyzing both global MHD simu- 
lations, and also simpler toy accretion models. The ultimate goal is to produce concrete 
predictions that can be compared directly with the large and continually growing body of 
X-ray observations of accreting black holes. 

Precisely because of the large interest in this topic, we give here a comprehensive de- 
scription of the technical components of our radiation transport code Pandurata. We hope 
that the techniques outlined below will be valuable to others who are interested in developing 
similar (or even better, more powerful) tools. We also present the results from a suite of 
simple numerical tests to verify the code, thus lending support and increasing our confidence 
in earlier work based on Pandurata. 



2. LOCAL ORTHONORMAL FRAMES 

The most general input for Pandurata is a body of tabulated data including the ex- 
trinsic fluid variables density, temperature, magnetic field, and the 4-velocity at each point 
in a three-dimensional volume. Multiple data slices in the time coordinate allow for studies 
of variability. The coordinates are Boyer-Lindquist for a Kerr BH with mass M and dimen- 
sionless spin parameter a/M . The fluid variables are given in physical cgs units for a specific 
BH mass and accretion rate. The source of the data is quite general, and Pandurata has 
been used successfully to analyse simulation data from the relativistic MHD codes HarmSD 



-4- 



dNoble et allboill : ISchnittman et al.ll2012l) and GRMHD JSchnittman et al.lboOGf ). as well as 



2D hydro simula tions (ISchnittman fc Rezzollal 120061) and analytic models for the accretion 



disk and corona (ISchnittman &: Krolik 



2009 



20101 ) 



We adopt a ( — h ++) metric signature, and convention where Greek indices run from 
to 3, and Latin indices are restricte d to spatial coordinates 1 to 3. The coordinate metric 
is given by ( iBoyer &: LindquistI 119671 ): 
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This allows for a relatively simple form for the inverse metric: 
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2.1. Simulation Data 



We briefly describe here the format of data from the HarmSD MHD simulations. Similar 
data can be generated from GRMHD, and any analytic model can be u nderstood as a subs et of 
the full tabulated simula tion data. As described in greater detail in iNoble et al.l ( 2011 ) and 
Schnittman et al.l (|2012[ ). the first step in post-processing the simulation data is to convert 
from code units of density and local dissipation to cgs units of density and temperature. 
Given the density everywhere, we integrate the optical depth along paths of constant (r, 0) 
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coordinates starting from both ^ = and 6 = n to get rtop{r, 6, (j)) and T}^ot{r, 0, (p). The disk 
midplane can then be defined as the surface 6'mid(r, 0) where rtop(r, ^^mid, 4>) = Thoti^, Omid, 4>)- 
When r(r, 6'mid5 <t>) > 1; the disk is optically thick and we define a top and bottom photosphere 
0(r, 0) such that rtop(r, 0top, 0) = Tbotl?^, ©bot, = 1- In Fi gure [J we show a slice in the 
(r, z) plane of simulation data from the HarmSD "ThinHR" run (INoble et al.ll2010l ). The local 
temperature is represented by the logarithmic color scale, and the contours show surfaces 
of constant r. In Figure [2] we show a three-dimensional representation of the photosphere 



surface 6 



top I 



for the same simulation data. 



From the photosphere surfaces, thermal photon are launched into the optically thin 
corona above and below the disk. Because the opacity within the disk is usually dominated 
by electron scattering, the seed photons are emit t ed wi th the limb-darkening and polarization 
dependence on angle given by IChandrasekharl (Il960l ). The spectrum is that of a diluted 
blackbody with temperature T^s and hardening factor /: 



with the black-body function. We take / 
effective temperature is given by 

Teff(r,0) = 



(4) 



1.8 (IShimura fc Takaharalll995l ) and the local 



1/4 



(5) 



where ITir^ (p) is the total integrated luminosity in the optically thick part of the disk (the 
factor of 2 comes from the fact that the flux is emitted equally from the top and bottom 
photospheres). As shown in Figure [H the gas has a constant temperature inside the disk for 
a given (r, 0), due to the high level of thermalization caused by the large optical depth. 

Synchrotron and bremsstrahlung seed photons can also be generated in the coronal re- 
gions, in which case we use an unpolarized, isotropic distribution function for the emission 
angles, as measured in the local fluid frame. Due to the high temperatures and low den- 
sities of the coronal regions, the net power in the coronal seed photons is typically much 
l ower than that of invers e- Compton scattering from the thermal seeds coming from the disk 
( jSchnittman et al.ll2012l ). 



2.2. Photosphere tetrads 



We begin with a short discussion of notation. As stressed in lMisner et al.l (jl973[ ). vectors 
are invariant geometric objects independent of coordinate system, and we represent them 
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Fig. 1. — Fluid density {top) and temperature (bott om) profile for a slice o f HarmSD data in 
the (r, z) plane, taken from the ThinHR simulation ( iNoble et al.ll2010l . l201ll ). Contours show 
surfaces of constant optical depth with r = 0.01,0.1, 1.0. Fiducial values for the black hole 
mass M = lOM© and luminosity L = O.lLEdd were used, as described in ISchnittman et al. 
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with bold font u, while the components in a specific basis are represented with italics: u^. 
We adopt a naming convention such that the components of a vector in the coordinate basis 
are represented by fi and in the local fluid frame by fi. The basis vectors themselves are 
labeled with (yu). For example, the coordinate basis is spanned by the vectors e(^) with 
components Cq^-, = ^q^-, , where 6 is the usual Kronneker delta. Note that the coordinate basis 
vectors are not normalized, and not even orthogonal in the Kerr metric: 

e{^) ■ e(^) = 9apef^)e1^^ = 9^u ■ (6) 



Einstein's Equivalence Principle, one of the bedrocks of general relativity, states that 
an orthonomal basis (a "tetrad") can be defined at any point in space. In fact, an arbitrary 
number of tetrads can be defined at any point, and are all related by Lorentz boosts and/or 
rotations. One particularly useful tetrad iri the K err metric is that of the zero-angular- 
momentum observer (ZAMO; iBardeen et al.l ( 119721 )). We denote the ZAMO frame with jl 
labels. It can be constructed from the coordinate basis by: 
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Any vector can be represented by its components in different bases: 



u 



and the components are related by a linear transformation E ~: 
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In our example of the ZAMO frame, E^- is given by 
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At each point on the photosphere we define a tetrad in the comoving fiuid frame (desig- 
nated with sub/superscripts fi) such that the time coordinate is in the direction of the fiuid 
4- velocity: 

= (11) 
In our notation, this equation says that the 4- vector tangent to the world line of an observer 
moving with the fiuid can be expressed in the Boyer-Lindquist coordinate basis with com- 
ponents /i, or in the local frame with components fi with e^^^^ = [1, 0, 0, 0]. The spati al basis 
vectors in the fiuid frame e^,^ are constucted via a method similar to lBeckwith et al.l (120081 ) . 
including a slight modification to ensure the right-handedness of the basis such that 6(5) is 
in the —9 direction. For completeness, we reproduce those definitions here: 
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(12a) 
(12b) 

(12c) 

(12d) 
(12e) 

(12f) 
(12g) 

(13a) 
(13b) 
(13c) 



From this tetrad basis, any other tetrad in the fiuid frame can be constructed from a 
simple rotation of the spatial basis vectors. We take as our preferred basis (now labeled 
with e(^)) one in which e(g) is normal to the photosphere surface. Whether we are using 
simulation data or an analytic model for the disk surface, the photosphere is described by a 
two-dimensional surface on the top and bottom of the disk: Qtop{^,'P) and 6bot('^5 0)- From 
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these functions, at each point on the photosphere we can construct two vectors tangent to 
the disk surface through the following process. Start with the coordinate-based vectors 

dQ 

dr'' = \0, Ar, —Ar, 0] (14a) 
or 

and 

dr = [0,0,^A<p,A<p], (14b) 

where Ar and A0 are the differential sizes of the fluid cell in questioij^. Next, subtract off 
the components parallel to e^^^: 

dr' = dr - (dr • e(£))e(£) (15a) 

and 

d0' = d(/) - d0 ■ e(£))e(£) . (15b) 

When dr' and dcp' are projected into the fluid frame, they will have only spatial components 
and will be tangent to the photosphere. In this basis, we can easily construct the normal 
vector by taking the 3- vector cross product: 

dz'' = e\.dr''d/^ . (16) 

This procedure has the added advantage of giving the proper area of the photosphere patch 
subtended by the vectors dr' and d(f)' by dA = |dz|. This formula for dA will be helpful for 
determining the amplitude of emitted flux from each patch of the disk, since the emission 
function is typically deflned in the local fluid frame. Because dr' and d0' are not generally 
orthogonal, we also deflne the dx and dy tangent vectors by dx = dr', dy'^ = and 

dy'' = A.d2^dx^ . (17) 

To complete the tetrad, we simply need to normalize the differential basis vectors. 
Returning to the coordinate basis, we have: 

ef^ = e^^^=u^ (18a) 

ef_^ = dx>'/{g^^dx''dx^f'^ (18b) 



er.^ = ±dy^/{g^pdy''dy^Y'^ (18c) 



e7.^ = ±dz^'/{go,^dz''dz^f'^ . (18d) 



The ± in the deflnitions for e(y) and 6(2) are chosen for the top (+) and bottom (— ) photo- 
sphere surfaces so that the spatial basis vectors are oriented in a right-hand fashion and to 
ensure that 6(2) points away from the disk body. In Figure [2] we show how these tetrad basis 
vectors are oriented on the photosphere surface. 



^For example, the ThinHR simulation uses Ar/r — 0.004 and A0 — 7r/128. 
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Fig. 2. — Three-dimensional representation of the disk photosphere surface Gtopl?", 0), along 
with the local spatial tetrad definitions of equation fllSap . The simulation data is the same as 
in Figure [H The color scale is a linear representation of the disk's local thermal temperature. 
The labels (r^, (pj) correspond to coordinates of the computational grid boundaries. 



In addition to launching photons from the disk surface, we often want the option of 
including seeds from within the corona, due to thermal bremsstrahlung, cyclo-synchrotron, 
or other radiation processes. Analogous with the comoving surface element defined above 
for disk emission, for coronal emission we need to define a volume element and associated 
tetrad at each point in the simulation space. Like the tetrads defined above, the time axis is 
defined along the local fluid 4-velocity u^- However, unlike the surface tetrads, the volume 
tetrads have no preferred orientatioqj, so we can simply use the spatial coordinate vectors 
projected onto the space orthogonal to u^: 




2.3. 



Coronal tetrads 



[0,Ar,0,0], 
[0,0,A^,0], 
[O,O,O,A0], 



(19a) 
(19b) 




•^For some specialized emission models, such as optically thin synchrotron, it may be convenient to choose 
a special orientation, e.g., with the e(^z) basis rotated to lie along the local magnetic field vector. 
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dr' = dr - (dr ■ e(£)) , (20a) 
de' = d6> - (d6> • e(£)) , (20b) 
d0' = d</> - (d</) • e(£)) . (20c) 

(20d) 

The proper volume element subtended by these vectors is given by the 3- vector triple product 
in the local fluid frame. While there is no real preferred orientation for the spatial axes, we 
still need to go through the process of defining some orthonormal basis to project the above 
vectors and thereby calculate vector products. In practice, we set e(i) along the dr' direction: 

clxf" = dr'", (21a) 

then set the y-axis roughly along the coordinate direction 



dy^ = e\.dx 'de^ , (21b) 



and the z-axis normal to both: 



dz^ = e\.dx 'dy^ . (21c) 
As above for the photosphere tetrads, the final step is to normahze all the basis vectors: 

<t) = (22a) 

ef-) = dx''/{g^pdx''dx^Y''' (22b) 

ef-) = dy^/ig^pdy'^dy^f/^ (22c) 

ef,) = dz"/{g^pdz''dz^Y'\ (22d) 

Unlike the photosphere case, since there is no "top" or "bottom" in the corona, we need 
not be concerned about the orientation of the 6(2) vector, and simply require a right-handed 
(x,y,z) convention. 



3. RAY-TRACING 
3.1. Geodesies 



The ray-tracing portion of Pandurata integrates the geodesic trajectories of photons in 
the Kerr metric. From the tetrad frames defined in the previous section, the initial direction 
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of a photon is selected from an isotropic distribution in the emitting fluid frame (limited to 
a hemisphere in the case of an optically thick photosphere surface). 



The geodesic integrater is the same as that described in ISchnittman fc Bertschinger 



(120041 ). based on a Hamiltonian approach. Because the Kerr metric is stationary, the mo- 
mentum conjugate to the time coordinate t is conserved, and corresponds to the (negative) 
specific energy of a particle (m^ = 1) or photon (m^ = 0). We can replace the affine 
parameter with the coordinate time and write the Hamiltonian as 
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with equations of motion 





dHi 


Itt ~ 


dpi 


dpi 


dHi 


dt 





(23) 

(24a) 
(24b) 



In Boyer-Lindquist coordinates, the Hamiltonian can be written thus: 



H{r,e,^,Pr,Pe., 



ujp^ + a{ ^pI + \pI + -^pI + rr? 



1/2 



(25) 



using the same notation defined above in equations fl3aH3el) . Because the metric, and thus 
Hamiltonian, is axisymmetric, Prp is also an integral of the motion. We are thus left with 
five coupled first-order ordin ary differentia l equations for {r,6,(f),pr,pe)- The third integral 
of motion. Carter's constant (jCarterlll968l ) 



Q=pI + cos^ e [a'^{m^ - pI) + pj/ sin^ 6] , 
is used as an independent check of the accuracy of the numerical integration. 



(26) 



For the numeri cal integration of geodesies, we use a 5*^-order Cash-Karp algorithm with 
adaptive step size (jPress et al.lll992l ). In Figure |3] we show the accuracy of the integrator 
by plotting the average deviation in the Carter constant for a selection of photons around a 
black hole with a/M = 0.99, as a function of step segments. We typically set the tolerance 
at 10~* per step, which we find allows sufficient sampling of the fluid variables near the black 
hole. Because of the frequent table look-ups required when using simulation data, there is 
little to be gained by using more advanced integratio n techn iqu es such as Bulir s ch-St oer or 
the semi-analytic approaches of iRauch &: BlandfordI (Il994j ) or iDexter &: Agoll (|2009[ ) that 
calculate the geodesic endpoint in a single integral evaluation, and are more appropriate for 
vacuum transport. 
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Fig. 3. — Convergence of the geodesic integrator, as determined by the accuracy of the con- 
served quantity Q. As expected, we find S^'^-order convergence for the Cash-Karp adaptive 
time step integrator. 
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N steps 



1000 



3.2. Polarization 

Pandurata is also capable of polarized transport along geodesies. The polarization 
vector is a space-like vector normal to the photon direction. For a photon with wave vector 
k, the polarization vector f is constrained by f • f = 1 and f ■ k = (jConnors et al.lllQSd ). The 
vector f is parallel transported along the geodesic path: Vkf = 0, but instead of explicitly 
solving the parallel tra nsport equation, we can take advantage of the com plex- valued Walker- 



Penrose constant K^p (IWalker fc Penrose! 1 19 70l : IConnors fc Starklll977l ) 



After solving for the wavevector k'^ along the geodesic path, K^p is given at any point 



by 



wp 



[{k'f - k^'f) + a sin^ eik^'f - k^f) 

-i[(r2 + a^){k'^f - f^k^) - a{k'f - k^ f)] sinO] (r - mcosi 



(27) 



Combined with the normalization factors f -f = 1 and f -k = 0, we have four linear equations 
for the four components of f'^. Because k is a null vector, we can always redefine f by a 
multiple of k: f ' = f -|- Ak, and thus write the polarization vector as 



f = [0, cos ipe\ + sin ipel] 



(28) 
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for some space-like basis vectors ei and e2 normal to k. 

The degree of polarization 5 < 1 is invariant along the ray path. When interacting with 
a distant detector or scattering off an electron in the fluid frame, it is convenient to employ 
the classical Stokes parameters /, Q, and U. In the (ei, 62) basis, we can write 

X = Q/I = Scos2ip, (29a) 
Y = U/I = 6sm2ip. (29b) 

One of the main advantages of this approach is that the Stokes parameters for each photon 
can simply be added at the detector, quite useful in a Monte Carlo calculation. Furthermore, 
/(z^), Q{v)i and U{v) can all be written in units of spectral density, which is the standard 
observable for many real detectors. 

For photons e mitted at an ang l e ^em to the normal of a scattering-dominated surface. 



we use the results of I Chandr asekharl ( 1l960l ) to get the initial polarization amplitude [ranging 



from ^(6*6111 = 0°) = up to 5{6em = 90°) = 0.12] and direction (parallel to the disk surface). 



3.3. Photon packets 

Because the geodesic photon trajectories are independent of photon energy, we can sig- 
nificantly improve the efficiency of the Monte Carlo calculation by tracking large numbers 
of photons simultaneously, covering a range of energies. We call the se computational ent i- 



ties "photon packets," which are analogous to the "superphotons" of iDolence et al.l ( l2009l ). 
except for the fact that theirs are monoenergetic and ours are broad-band. We also assign a 
single polarization angle and degree to the entire photon packet. This is an approximation 
that works well for vaccuum transport and coherent scattering, but will break down when 
including scattering at high energies hv > me(? as the electron cross section becomes more 
energy- dep endent . 

Each photon packet is weighted by a number of geometric emission factors. For example, 
a photon packet emitted from a small patch of optically thick, scattering-dominated accretion 
disk would have a spectrum of 

FT = ^B^if Tes) /limb(^em) COS ^em dA , (30) 

where is a function that has units of spectral luminosity [erg/s/Hz]. Here / is the same 
hardening function introduced above in equation (jlj), cos 6'em is a geometric fa ctor for emission 



from an optically thick surface, /umb is a limb- darkening function given by I Chandr asekhar 
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( ll960l ). dA is the proper area of the emission region [see eqn. ( !T6|) above], and dQ = 2ti /N^\i 
is the proper sohd angle of a hemisphere sampled evenly by A^ph photon packets. Lastly, 
1 / = dr/dt is a. relativistic correction factor to convert from time in the emission frame to 
that in the coordinate or distant observer frame. 

In order to account for the spectral redshift, we store both F and z/ at a set of discrete 
points. When transforming from the emitter to observer frames, F is invariant (units of s~^ 
and Hz~^ cancel)0, while u transforms as follows. If the photon packet is emitted in a frame 
with fluid 4- velocity ■u'^(em) and photon 4-momentum fc^(em), and observed in a frame with 
M'^(obs) and k^{ohs), then we can write the redshifted frequencies as 

M^(obs)/c^(obs) 

z/(obs) = z/(em) — — - — - . (31) 

Whenever the photon packet scatters off the disk or an electron in the corona, the frequencies 
Vi are updated and the old "observed" frame becomes the new "emitted" frame. When 
the photon packet reaches an observer at infinity, M'^(obs) = [1,0,0,0] and the well-known 
redshift relation is reproduced. 

For this distant observer, the angle of polarization tl) is measured by projecting f onto the 
(61,62) = (e<^, —eg) basis. For an observer oriented with the black ho le spin axis projected 



in th e verti cal direction, ^ = corresponds to horizontal polarization ( ISchnittman fc Krolik 



20091 . 12OIOI ). Given ^, 5, and Fy^ the spectral luminosity form of the Stokes parameters are 
simply 

Qy = FJcos2il), (32a) 
ily = FJsm2ip, (32b) 

where Q and U are related to the Stokes parameters Q and U hj a factor of (F/I). After 
summing over a large number of photon packets, we then invert equation (!32|) and return to 
the Si^u), ifj^u) representation. 



3.4. Emission and absorption 

Along each geodesic path, we can also include local emission and absorption pro- 
cesses such as bremsstrahlung or synchrotron. This is the predominant method for gen- 
erating light curves and spectra in codes that shoot photons backwards from a distant 



"'For a discrete function Fi, the number of photons emitted per coordinate-frame second between Vi and 
Ui + dvi is Fi dvi/ ihvi), where h is Planck's constant and Ui are measured in the local emission frame. Because 
Vi and dui transform the same under Lorentz transformations, Fi is invariant. 
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2006: Noble et al 



2007 



2009 



observer ( Broderick fc Blandford 2004 ■ Schnittman &: Bertschingerl |2004J : ISchnittman et al. 



Dexter fc Agol 



20091 ) ■ In the fluid frame, the radiation trans- 



port equation is given by 



ds 



(33) 



where ds is the differential path length and I^, j^, and a,^ are respectively the spectral 
intensity, emissivity, and absorption coefficient of the local ffuid. The absorption coefficient 
is related to the opacity k^, through the density p: = pn^. Defining the optical depth Ti, 
through 

dT,y = a^ds, (34) 

the transfer equation can be written as 



— Oy 1, 

dTy 

where the source function is defined as Sy = ju/dy. 



(35) 



Both ly and 5*;^ have the same properties under Lorentz transformations, namely ly/u^ 
and Sy/u^ are both invarian t. Other invariant terms are the optical depth Ty, uay, and jy/i'^ 
( iRybicki fc LightmanI 120041 ) . Thus if we can solve the non-relativistic radiative transfer 
equation ( l33l) in the local fiuid frame, then in any other inertial frame (e.g., the ZAMO 
tetrad), the special relativistic version can be written 



dly f^V ; \ > T 



(36) 



Here the fluid frame (where jy and ay are deflned) is the primed frame, and the "lab frame" 
unprimed. 

The above analysis, while quite useful for special relativistic flows in the locally flat 
ZAMO basis, ignores all general relativistic effects of curved spacetime around the black 
hole. To include these effects, we need only shift the frequencies Vi from one geodesic step 
to the next, due solely to the gravitational redshift, and we can treat each computational 
step as a new observer relative to the previous step, as in equation (j3T|) . 



4. SCATTERING 

We allow for two types of scattering in Pandurata: Compton scattering off free electrons 
in the corona, and scattering off an optically thick disk (which in turn is characterized by 
repeated scatterings in the relatively cool atmosphere). Because electron conserves photon 
number, our photon packet approach is ideal for modeling these processes. 
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4.1. Coronal Scattering 

The first step in the scattering process is to determine whether a scattering event takes 
place at all. To do this, we transform into a local inertial "lab" frame, generally taken to be 
the ZAMO frame discussed above in Section [2?2l In this frame, the photon moves a distance 
of dp = rfi'-dx^dx^ in a single geodesic integration step dt. Then the total optical depth to 



scattering along the path is 



dr = dl Kes Plab = dl Kcs Pfluid^^^ , (37) 

^'lab 



where the last equality comes from the invariance of z/ctj, (IRybicki &: Lightmarul2004j ). with 
the absoption coefficient a^, = n^sP for electron scattering opacity. Given dr (typically much 
less than unity), the probability of scattering is P = 1 — e^*^"^. 

When a photon does scatter off a free electron, we carry out the scattering calculation in 
the electron's rest frame. This requires two coordinate transformations: from the coordinate 
basis (denoted with /i super/subscripts) to a fluid-frame tetrad (/i), and then a Lorentz boost 
from the fluid frame to the electron's rest frame (/2'). The transformation from coordinate 
basis to corona fluid frame is the same as given above in Section 12. 3[ In the fluid frame, the 
electron velocity /3 = f/c is taken from an isotropic Maxwell- Juttner distribution 



where 7 = 1/ yl— 6t = kT/meC^, and K2 is the modified Bessel function. See Appendix 
[B] for a description of our algorithm for generating a Monte Carlo sample of velocities that 
satisfy equation (l38l) . 



Following iMisner et al.l (119731 ). we construct a generic Lorentz boost in the direction of 
the electron 4- velocity u'^ = ['y,'yPn^]: 

= [7,7M (1^1 = 1), 

= 7, 

A*5 = A^; = -/37n^ 

A^'^ = A^' = (7 - l)n^n^ + 5^~K (39) 

The photon momentum in the electron frame is thus given by p'^' = A^^p^. 

Without loss of generality, we can carry out one more transformation and define the 
initial photon direction to lie along the z-axis in the electron frame. The x-y plane is 
decomposed into Si and £2-, where the initial polarization is aligned with Si. The scattered 
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Fig. 4. — Schematic of the scattering geometry in the electron frame. The incoming radiation 
is polarized along the e\ direction. The scattered radiation kj makes an angle G with e\ 
and Q with kj. When projected onto the e\ — £2 plane, k/ makes an angle with £\. 




radiation kj makes an angle B with e\ and Q with kj, as shown in Figure HI For unpolarized 
incident light, we can define e\ to lie in the plane of kj and kj, with O + = 90°. 

For photons polarize d along e^ , the angle-deperi dent cross section aiff) is given by the 
dipole scattering formula ( iRybicki fc Lightmanll2004j ): 

da 



Tq sin^ 6 



Tq cos^ 9 cos^ I 



(40) 



pol 



where (p is the standard azimuthal angle measured with respect to £1. Here the classical 
electron radius tq is given by 

^2 



2.82 X 10"" cm. 



(41) 



For photons scattering in the kj-£2 plane, the cross section is constant: da^Q = 7T/2)/dQ = 
Tq. For unpolarized light, we define £1 as lying in the scattering plane, so the scattering 
angle with respect to £2 is tt/2. Because unpolarized light is an equal combination of £1- and 
£2-polarized photons, we can reproduce the familiar cross section for unpolarized scattering: 



da 

dn 



unpol 



(daim 

\ dn J 



1 

2 

1 2 / -1 2 

-ro(l + COS 



+ 



pol 



da{n/2) 

dn 



pol 



(42) 
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Fig. 5. — Definitions of polarization axes in pre- and post-scattering coordinates, kj, kj, sy, 
and £m are all in the same plane, while £± and e'^ are normal to that plane. 




For an arbitrary polarization degree 5, the cross section can be written as the sum of 
unpolarized light with weight (1 — 5) and purely polarized light with weight 5: 

^ = ^r2(l-5)(l + cos2^) + r25(l-sin2^cos2 0) 

= ir2[(l + cos2^) -(5sin2^cos20]. (43) 

Given the angle- dependent cross section, we can either pick the scattering angles (6', 0) 
directly from a distribution function derived from ( 143 p . or alternatively, pick the angles from 
a uniform distribution, and give the scattered flux a weight based on the cross section. We 
compare these two methods in Appendix [Cl 

Once the new photon direction is detemined, we need to calculate the angle and degree 
of the post-scat t ered p olarization. Here we follow the Rayleigh matrix method described in 



Chandrasekhad (jl96d ). We define yet another coordinate system with £3 parallel to kj, e\\ 
in the scattering plane defined by kj and kj, and £_l normal to that plane. Likewise, we 
define a post-scatter frame with £3 parallel to kj, £'j_ = £_|_, and Sy in the scattering plane, 
but normal to kj (see Fig. [S]). In this frame, the initial polarization vector can be written 
fj = cosipe\\ + smip£_^_ and the final polarization is fj = cosip'e'^^ + smip'e'^. 

The standard Stokes parameters are given by the intensity I, Q = 5/cos2?/', U = 
(5/ sin 2?/', and V = (electron scattering never leads to circularly polarized light). We 
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further define 



^ i(/ + g) = i(l-5)/ + 5/cos2^ 
= 1{I-Q) = ^-{1-6)1 + 6Ism'^ 



I = [Ii\,IuU,V] 
and the Rayleigh scattering phase matrix 



R 



(44a) 

(44b) 

(44c) 











\ 





1 














cos^^ 





V 








COs6' y 



Then the scattered Stokes parameters are given simply by I' = RI, /' = 
Q' = /|| — Note that the cross section fHSl) can be reproduced by writing 

1 , 



(45) 



+ Ji, and 



I' = cos^ eh + 



5)1(1 + cos^ 6) + 5I{cos^ 9 cos^ ip + sin^ ip) 



givmg 



/' 1 

_ = -(1 -5)(1 + cos2^) + 5(1 -sin^^cos^^) 



(46) 
(47) 



now with ip taking the place of from equation 
Lastly, f/ is constructed by 



6' 
f/ 



VQ'^ + 
/' 

^tan-i(t/',Q') 
cosV^'^N + sin 



(48a) 

(48b) 
(48c) 



At this point, the polarization vector and photon 4-momentum are transformed back into 
corona fiuid frame, then to the coordinate frame, and then the geodesic propagation continues 
as before, until the photon packet scatters again, is absorbed by the black hole, or reaches a 
distant observer. 

During this scattering process, the photon packet's array of frequencies had to be ad- 
justed three times: once when transforming from the fiuid frame to the electron rest frame, 
once when losing energy to the electron recoil, and once when transforming back to the fiuid 
frame. The first and last transformations are simple Lorentz boosts, and the frequency scales 
like the photon energy: u'/u = p^' /p'^, with p*' = A*'^'^. For the scattering losses, we need 
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to scale the frequency bins such that the number of photons in each bin is conserved, while 
losing energy according to the Compton recoil formula: 



E 



f 



Thus the frequency scales like 



V 



l + -^^(l-cos^) 



hv , 

1 H ^(1 - cos I 

mpC^ 



(49) 



and the size of each bin scales like 

dv' _ 
dv 



hv , 
1 H ^ 1 - cos I 



-| -2 



(50) 



(51) 



The number of photons in each bin dN^, = F,jdv/{hv) is conserved in the scattering event, 
so we find that the effect of Compton recoil on the spectral luminosity is 



F' 

V 



1 + 



hu 



[1 — cos I 



(52) 



At very high energies hu 3> rUeC^, this leads to a "pile up" of photons and large peaks in 
the photon packet spectrum. In reality, this effect would be mitigated by incorporating Kline- 
Nishina cross sections, which decrease with energy, yet are incompatible with our photon 
packet approach that treats all photons as identical regardless of frequency. In practice, we 
are generally interested in problems where the characteristic photon energies are significantly 
below rrieC^, so the photon pile up is rarely an issue. 

While some energy is lost to Compton recoil in the electron frame, the more typical 
effect is inverse- Compton scattering, where energy is transferred from the electrons to the 
photons. For electrons with Lorentz factors 7 in the fluid frame and low-energy photons 
with hu/rrieC^ ^7^ — 1, the ratio of energies of the photons b efore scattering, in the res t 



frame of the electron, and after scattering is roughly 1:7:7^ (IRybicki fc LightmanI 12004 ). 



For coronal electrons with temperature ~ 140 keV, low-energy seeds will, on average, double 
their energy after every scattering event, making inverse- Compton a very efficient radiative 
process. 



4.2. Disk Scattering 



At each step along the geodesic trajectory, we determine whether or not the photon 
packet has crossed the photosphere surfaces 6top('", 0) or 0bot(^, 0)- If it has crossed this 



- 22 - 



boundary, we follow a procedure similar to that described above for coronal scattering. First, 
we use the conserved quantities /t^p, f ■ f = 1, and f ■ k = to solve for the polarization 
vector f in the coordinate frame. Then we transform f and k into the local fluid frame of 
the photosphere tetrad e(^), with e(5) normal to the disk surface, as in equation f llSap . 

In this frame, the scattering off the disk surface is calculated using the analytic expres- 
sions for refle ction off a diffus e semi -infinite plane, derived by Chandrasekhar and given in 



table XXV of IChandrasekharl (jl96d ). As in equation fl44p above, we can write the incom- 
ing photon beam as a vector of Stokes parameters for the flux I = [/||,J^,f/] {V = for 
linearly-polarized light, the only relevant case for our scattering-dominated systems). Then 
the outgoing intensity is given by 

//^\ 

I'(/i,¥') = 



\u'i 



4/i/io 




QS(/i, /io, ¥?o) MX , (53) 



where (/zq, V'o) are the incident angles in the fluid frame with /iq = \kz\, (/i, ^) are the outgoing 



angles , and Q and S are the transfer matrices defined in Section 70.3 of I Chandrasekhar 



(Il960l ). Unlike the coronal scattering case, where we use the differential cross section P3|) 
to determine the post-scatter angles, for diffuse reflection off the disk, we simply choose a 
random angle (/i, if) from a uniform distribution and then weight the outgoing intensity by 
/'// from equation (l53l) . Thus any individual reflection does not conserve photon number, 
but the angle-averaged process does. From I', we are able to reproduce 5', ijj' , and thus f 
and k' as above, which are then transformed back into the coordinate frame and continue 
their geodesic propagation through the corona. 

This method for diffuse reflection can be checked against coronal scattering experiments 
where we scatter incoming photons off a semi-infinite plane of free electrons. We find excellent 
agreement between the analytic and numerical approaches, as shown below in Section [51 

As with the coronal scattering, high energy photons can lose energy to Compton recoil off 
the electrons in the cool disk, leading to the reflection hump seen in many AGN observations. 
While this process is technically angle- dependent, as a simplification, we average over all 
incoming and outgoing angles, as well as the number of individual scatterings typically 
responsible for diffuse reflection {N^c-At ~ 3), and use the recoil formula 

p' f hu 

-=1 + 3--^ • (54) 



This energy lost by the photons can then be reprocessed by the disk and emitted at thermal 
energies. 
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Fig. 6. — Comparison of observer-to-emitter and emitter-to-observer ray-tracing paradigms 
for a relativistically broadened emission line. The disk extends from i^out = 15M all the way 
to the horizon, and the emissivity profile scales like r~^. 
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NUMERICAL TESTS 



Here we present a number of test problems to verify Pandurata's accuracy and reliability. 
We begin with vacuum transport of geodesies from the disk to a distant observer. To test the 
tetrad construction methods outlined in Section \2.2\ we calculate the relativistic broadening 
of iron lines from a thin disk around a Kerr black hole, comparing the emitter-to-observer 
and obser ver-to-emitter paradigms. The observer-to-emitter approach is well-known in the 
literature ( iRauch fc Blandfordlll994i ISroderick fc Blandfordll2003l . 12004 ). It is also relatively 
straight-forward conceptually, since it doesn't require the use of any tetrads or proper area 
calculations. One simply shoots rays backwards from a distant observer, and integrate the 
geodesic path unt il the ray crosses the midplane, where the fluid 4- velocity can be determined 
analytically as in iNovikov fc Thornd ( 1l973l ). This gives the redshift of the emission line as 
seen by the observer, and the spectrum is given by the invariant ly/v^. 

In Figure O we show the shape of a relativistically broadened emission line as viewed by 
observers at different inclination angles for the spin values a/M = and 0.99. In all cases. 



the emissivity profile scales like / 



and the outer disk is truncated at r = 15M. The 



disk extends all the way into the horizon, with the fluid velocity inside the ISCO determined 
by conserving the energy and angular momentum at the ISCO, and solving for u"^ from the 
relation m^m'^ = — 1. For the observer-to-em i tter c alculation, we use the same ray-tracing 
code described in ISchnittman fc Bertschingerl (120041 ). with 10^ photons evenly spaced in the 
image plane for each inclination. We find excellent agreement in all cases, validating our 



Fig. 7. — Error estimate e as a function of photon number for a relativistically broadened 
iron line, defined in equation f l55|) . 40 inclination bins were used. 




emitter-to-observer techniques, at least for planar test-particle orbits. 



This test in turn naturally leads to a simple convergence test for our Monte Carlo code. 
Integrating over energy and observer inclination angle i, we can apply the following metric 
to estimate the error due to the use of a finite number of photons: 



[JdcostJdE{lUE)-lUE)y] 



1/2 



[J d cos^ J dEPjE)] 



1/2 



(55) 



with I\o{E) the spectrum calculated at low resolution, compared to the theoretically perfect 
spectrum Ihi{E) calculated at high resolution. As expected for a Monte Carlo calculation, 
we find that the total error scales with photon number like iV~^/^, as shown in Figure [71 
This is consiste nt with similar spec tral calculations done with the Monte Carlo radiation 
code grmonty ( iDolence et al.l 120091 ). Also shown in Figure [7] are the errors ^(A^) for the 
observer-to-emitter approach, using a total of 40 inclinations for both cases. We note that 
the emitter-to-observer method is more than a factor of two more efficient for the same 
calculation. This is because we can selectively shoot more photons from the inner regions of 
the disk, but in the reverse method, the photons are di stributed evenly in the image plane 
(this uniform distribution is not strictly necessary; e.g.. iNoble et al.l ( 120071 ) use an adaptive 
grid to improve resolution in bothros). 

The next test is similar, but also includes polarization effects. Instead of an emission 
line with J(r) ~ r~^, we use the diluted thermal spectrum for a Novikov-Thorne (NT) 
disk with an inner edge at the ISCO. The emission has the polarization and limb darkening 
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Fig. 8. — Comparison of observer-to-emitter and emitter-to-observer polarized images for a 
NT disk with polarization given by a scattering-dominated atmosphere. The disk extends 
from an inner edge at the ISCO out to i?out = 15M. The black hole has spin a/M = 0.99 
and the observer is at an inclination of z = 75°. The intensity color scale is logarithmic 
and the polarization vectors are linearly proportional to the local degree of polarization, as 
observed at infinity. 
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appropriate for a scattering-dominated atmosphere (jChandrasekharlll960l ) . For the observer- 
to-emitter approach, in addition to utilizing the l^jv^ invariance, we also parallel-transport 
two polarization basis vectors corresponding to the two axes in the observer plane normal to 
the photon propagation direction. Then, when the ray intersects with the disk, we calculate a 
local tetrad in order to determine the local angle of incidence and thus degree of polarization. 
The direction of polarization is projected onto the parallel-transported basis vectors to give 
the observed angle at infinity. 

The two approaches give identical results, as shown in the images in Figure [HI for a 
Kerr black hole with spin ajM = 0.99, i?out = 15M, and observer inclination angle 75°. 
The color code is logarithmic in total intensity and covers four orders of magnitude, and 
the small vectors scale linearly with degree of polarization. For the purposes of comparison, 
we have not includ ed returning radiation here, despite the import ant effect it has on the 
polarization signal (lAgol fc Krolikll2000l : ISchnittman fc KrolikI l2009l ) . In fact, it is precisely 
due to the critical importance of returning radiation that we w ere forced to employ the 
emitter-to-observer approach in ISchnittman fc KrolikI (120091 . l2010l ) . 

In Figure M we show the observables of polarization degree and angle as a function 
of energy for a range of inclination angles, assuming an Eddington-scaled accretion rate of 
rh = 0.1 and black hole mass M = IOMq. Again, we find excellent agreement between the 
emitter-to-observer and observer-to-emitter methods. 
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Fig. 9. — Polarization degree and angle as a function of photon energy for a black hole with 
spin a/M = 0.99, luminosity O.lLEdd, and mass lOM©. The disk extends from the ISCO 
out to r = 15M. The observer-to-emitter (solid curves) and emitter-to-observer (diamonds) 
methods agree closely over a range of inclinations. 




Next, we move on to testing the coronal scattering algorithms. We focus on a plane- 
parallel geometry with an optically thick disk covered by a corona with variable optical 
depth T and electron temperature Tg. In Figure [TO] we show the effects of a scattering 
atmosphere on the emergent flux and polarization as a function of angle. The seed photons 
are emitted isotropically from the disk surface with zero polarization, then scatter through 
a cold corona. Photons that scatter back to the disk are reflected via the diffuse scattering 
formula of equation (!53l) . In the l imit of r oo , we r eproduce the limb-darkening and 
horizontal polarization results from IChandrasekharl (119601 ). Table XXIV. 



In Figure [TT] we carry out a similar scattering experiment, but with the seed photons 
incident from above the disk along a sin gle direction. Setting r = 10, we should reproduce the 
analytic diffuse reflection expressions of I Chandr asekharl ( 119601 ). shown there i n Figures 24-25 
for an incident unpolarized beam with cos^o = 0.8, 0.5 and (fQ = 0. Followin g I Chandr asekhar 
(Il960l ). we plot the Stokes parameters /, Q, and U as a. function of reflection angle, normalized 
to the incident intensity. The asterisks are the Monte Carlo calculation, and the solid curves 
are the analytic predictions. 

On the left-hand side of each plot, we show the polarization as a function of 6 for 
ip — ipQ = 0°,±180°. The value of 9o is designated with a vertical dashed line. Negative 
values of 6 correspond to photons reflected back in the general direction of the incident 
photons, i.e., ip — (po = 0. Thus we see a natural peak in the intensity corresponding to 
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Fig. 10. — Outgoing intensity and degree of polarization for radiation emitted from a 
scattering-dominated atmosphere, as a function of inclination, and for a range of corona 
optical depths. The face-on orientation is unpolarized due to symmetry. 
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backscattering as in equation (l42l) . Similarly, the degree of polarization is maximized for 90° 
scattering, and oriented in the plane of the disk {Q > 0). 

On the right-hand side of each plot, we show the Stokes parameters for photons scattered 
with — ifQ = ±90°. In this case, the planar symmetry is broken and we find a non-zero 
value of U. Again, the degree of polarization is maximized for scattering angles near 90°. 

Lastly, we t est the inverse- Compt o n effe cts of a hot corona by reproducing the AGN- 
type spectra of iPoutanen fc SvenssonI (119961 ). The seed photons are again isotropic and 
unpolarized, with a blackbody spectrum with Tbb = 10 eV. When refl ecting off the cold disk, 
we irn plement the Compton recoil losses of equation ( 1541) . Following iPoutanen fc Svensson 
(119961 ). we also include atomic absorption in the disk with an extre n iely s imple toy model 
based on the photoelectric cross sections of Morrison &: McCammod (1l983l ). 



In Figure [12] the corona temperature is 56 keV with optical depth r = 0.5, and in Figure 
13) = 352 keV and r = 0.05, corresponding to Figures 5 and 6 in IPoutanen fc Svensson 
(jl996 ). In the upper panels we show the observed flux at two inclination angles cosi = 



0.11, 0.5, and in the bottom panels we show the polarization degree 6{%) = Q / 1 x 100 . In all 
panels, the solid curves correspond to the total flux, while the (dotted, dashed, dot-dashed, 
triple-dot-dashed, and long-dashed) curves represent subsets of the flux, binned by number 
of coronal scatterings (0, 1, 2, 3, > 5). Photon packets that return to the disk suffer pho- 
toelectric absorption and Compton recoil losses, and are then launched again from the disk, 
resetting Ngcat to zero. Thus the dotted curves in Figures [12] and [T3] have significan t powe r 
around the Compton hump at 10-100 keV. As discussed in ISchnittman &: KrolikI (120101 ). 
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Fig. 11. — Stokes parameters for scattering of unpolarized light off an op tically thick at- 
mosph ere. See text for description. Compare with Figures 24 and 25 of I Chandr asekhar 
fll960h . 
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more scatterings in a sandwich corona effectively constrain the geometry and increase the 
amplitude of polarization at high energies. 

We find excellent agreement overall, but are clearly dominated by Monte Carlo noise 
above ~ 100 keV. For these disk and coronal parameters, this corresponds to seed photons 
that have already scattered on average over 25 times, so it is very difficult to resolve any 
polarization signal at the few percent level. Additionally, due to our photon packet algorithm, 
we are limited to energy-independent electron cross sections, so we should expect that the 
accuracy of our spectral predictions breaks down much above 100 keV anyway. 



6. CONCLUSION 

We have presented the technical details behind the general relativistic radiation trans- 
port code Pandurata. Its capabilities include optically thin emission and absorption, Comp- 
ton scattering, polarization, spectral and timing analysis, and flexible geometries that allow 
analysis of numerous accretion models and MHD simulations. We have discussed a number 
of practical challenges that may also face other teams working to develop similar ray-tracing 
codes, such as the method of weights in the scattering kernel. 

This is by no means the final word on Pandurata. Its great strength lies in its flexibility, 
and we envisage numerous upgrades and improvements in the near future. These will include, 
but not be limited to, detailed ionization balance in the disk photosphere for improved AGN 
modeling, time interpolation between simulation snapshots for generating more accurate 
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Fig. 12. — Spectra and polarization of flux from a disk and corona with a slab geometry cor- 
responding to an AGN accretion disk, as described in the text. The solid curves correspond 
to the total flux, while the (dotted, dashed, dot-dashed, triple-dot-dashed, and long-dashed) 
curves correspond to A'scat = (0, 1, 2, 3, > 5). For cl arity, only the solid cu r ves ar e shown for 
the polarization degree. Compare with Figure 5 of iPoutanen fc Svenssoru (jl996[ ). 
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light curves, and the inclusion of more sophisticated emission and absorption processes (e.g., 
angle- dependent synchrotron) to model low-luminosity sources such as Sgr A*. Perhaps most 
important, we will work to close the final remaining gap between theory and observation by 
incorporating Pandurata spectra into a data analysis framework like xspec and making it 
publicly available to the X-ray astronomy community. 
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A. Hamiltonian Equations of Motion 

The equations of motion for tlie Hamiltonian H in Boyer-Lindquist coordinates, as given 
in Section [21 are repeated here for completeness: 



H{r, 9, 4),pr,pe,Pcf>) = -pt = ujp^ + a (^^p^ + ^pl + + m^^ 

and according to classical theory: 

dx^ dH 
dt dpi ' 

dpi dH 



1/2 
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For convenience of notation, we define tfie quantity D"^ as 



A 2 ,1 2 , 1 2 



m 



Then for an arbitrary variable y e the partial derivative of H can be written 

q;2 9L>2 



(A2) 



9i7 (9 "^"^/^ "'" 

dy dy dy 2pt + ujp4, dy 

The first set of Hamiltonian's equations are straightforward to produce: 



(A3) 



dr 
'di 

de 

It 
d^ 
Itt 



dpr 

dHi 
dpe 
dHi 
dp4> 



Pr 

Pt + a;p0 p2 
Pe o? 
Pt + bJP4, p^ 



u — 



P<p 



a 



Pt + cop^ 



(A4a) 
(A4b) 
(A4c) 



The momentum equations are a bit more involved, but there are only two of them (for p^ 
and p0] p^ is conserved): 



dpr duj Pt + up^ da a 
dt dr^'^ a dr 2{pt + ojp^) 
dpo doj Pt + uip^ da a^ 

~dt ^ ~ a^^^ ^ a 'm ^ 2{pt + ujp^) 



' d [A , 1 , 1 2 



(A5a) 
(A5b) 
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The relevant spatial derivatives are as follows: 



dr 
du! 

de 

da 

dr 
da 

ae 

dr 
d^ 

de 

d_ f J_ 
dr \ci7^ 
d_ fj_ 

m 

d_ / A 

Or 

d_ / A 

d_ ri_ 

dr 

d_ ( 1 

de 



2Ma 

,2 



2Ma 
1 da^ 

2a dr ' 
1 da^ 
2^~d9' 
'2M 



3r^ + 0^(1 + cos^ e)-^ cos^ e 



2Ma^ - aV - — 1 sin e cos e 
r 



— —a 



—a 



Ap2 



a — r 



^ 2r^a'^ sin^ e 



AMa^r sin 6^ cos e{a^ + r^) 



— J sin^ 6* ( r + 
4 sin ^ cos ^ 



Ap2 

2A#«- sill- e()s2 - ,-2) 



2 ^r-M- 
P \ P 



2Ma^ sin^ e 

rA 



+ 



+ 



+ (r^ + a') 



TO^Asin^cos^ , 

,4 ' 



2r 



ro sin e cos 6* . 



(A6a) 

(A6b) 

(A6c) 
(A6d) 
(A6e) 

(A6f) 
(A6g) 
(A6h) 
(A6i) 
(A6j) 
(A6k) 
(A61) 



B. Monte Carlo Sampling of McLxwell-Juttner Distribution 

For any normalized distribution function f{x) with x e (— oo, oo), one can always define 
the cumulative distribution function 



cdf(x) = F{x) = / f{x')dx', 

J —oo 



(Bl) 



with F{—oo) = and F{oo) = 1. Then by selecting a uniform random number A G [0, 1), 
the choice x = F~^{X) will be distributed according to f{x). However, in most cases, F{x) 
cannot be written in closed form, so other methods are required. 
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One simple technique described in iPress et al.l ( 119921 ) is the "rejection method," where 
an auxihary function g{x) is used, where g{x) > f{x) everywhere, and G{x) is easy to 
calculate. We begin by selecting a trial Xq = G'^^(Ao), then pick another random deviate 
Ai- If Ai < f{xo)/g{xo) then xq is selected as a representative sample of f{x), else we try 
again with a new Aq. Of course, if g{x) is large enough, it is easy to ensure that it is greater 
than f{x) everywhere. However, the efficiency of this method is limited by the ratio of the 
areas un der the two curve s f{x) and g{x), so it is desirable to pick g{x) as close to f{x) as 
possible (jPress et al.lll992l ). 



For the Maxwell- Juttner distribution defined in equation ([38 



fil) ~ 7^/3 exp( -7/6*7 



we choose an auxiliary function 



g(j) ~ 72exp(-7/6'' 



This gives 



G(7) 



1 



(B2) 



(B3) 



(B4) 



20T + 1 

for the cumulative distribution function. Inverting (]B4|) isn't trivial, but can be done nu- 
merically with a simple root finder. For these choices of /(7) and (7(7), we find excellent 
efficiency for this algorithm of ~ 90%. 



C. Comparison of Scattering Kernels 

As described in Section HI there are (at least) two different ways to implement the 
scattering of polarized light off of free electrons. 

The method of weights picks a random scattering angle from a uniform distribution of 
cos^ G [—1,1] and G [0,27r), then weights the scattered beam of photons by the cross 
section in that direction, normalized by the average cross section to conserve flux. By 
integrating equation (j^3|) over (j), this resembles the classical cross section for unpolarized 
light: 

uj{e) = j = ^{cos^e + i). (ci) 

Because repeated scatters tend to increase the level of polarization (indeed, in the mi- 
croscopic limit, every photon has 5 = 1), we will focus on the case where 6 = 1, giving 

u;(e,0) = ^(cos^^cos^^ + sin^^). (C2) 
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For angles uniformly distributed in cos 6 and (f), one can show that the probability distribution 
function (pdf) for w is 

1/ 2 



P{w) = -i^l--wj (C3) 

for < w < 3/2, and otherwise. 

For multiply-scattered photons, the weight function is multiplicative, since the individ- 
ual scattering events are uncorrelated. For n scatters, the net weight is given by 

n 

W = l[wi. (C4) 

i=l 

To determine the pdf P{W), we define a new variable Z: 

n n 

Z = InVT = ^Inwi = ^Zi. (C5) 

i=l 1=1 

For large values of n, the central limit theorem dictates that the distribution of Z should be 
Gaussian: 

1 exp(-i^-I^), (C6) 

where fiz and o"^ are respectively the mean and variance of P{z). From equation (1C3P and 
the variable transformation z = Inw, we have 

P{z) = P{w)^='4(l-le^y'^\ (C7) 



dz 3 \ 3 ^ 

with z G (— oo,ln3/2]. This gives /i^ = —0.208 and = 0.710. Now we see that the pdf 
P{W) is given by a log-normal distribution: 

P,W)^—L^e.JJ^^^^^^). (C8) 

For photons random-walking through an atmosphere of optical depth r, we find the pdf 
of the number of scatters required to escape can be closely approximated by 

Then the net distribution P{W) for all scatting orders is simply 



oo 



P{W]t)= dnP{n]T)P{W]n). (CIO) 
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The relative contribution to the spectrum from photons with a weight in the range 
(py, + dW) is P{W)WdW, so we require P{W) to decrease faster than W'"^ for large 
W if the calculation is to converge. In Figure [T3] we plot W'^P{W) for a range of r. Our 
analytic results suggest that for r > 2, any polarization spectrum formed using this Monte 
Carlo weighting method should be dominated by the rarest, highest-weight photon packets, 
confirming what we have seen in trial runs with large r. Now, in practice, the convergence 
is not quite as bad as Figure [H] suggests, for two primary reasons. First, the seed photon 
packets have little or no polarization, so the initial weighting function more closely resembles 
equation flCl|) . which leads to a significantly tighter range in W: yUunpoi = —0.027 and 
^unpoi — 0.047 (in this unpolarized limit, the weight method converges for all optical depths 
up to r > 200). Second, for the small-to- moderate optical depths of r < 5, the typical 
number of scatters n is still small enough that the mean value theorem does not strictly 
apply, in effect cutting off the high-weight tails in flC6p and further reducing the contribution 
from statistical outliers. 

However, for r > 10, the polarization of a typical photon bundle reaches 5 — )■ 1 after 
just a few scatters, and the large number of total scattering events allows us to reproduce 
these analytic results with numerical tests of the Monte Carlo code. In Figure [TBI we show 
the distribution of weights from a calculation using unpolarized seed photons, scattering 
through optical depths of r = 2, 5, 10. While we find that the r = 5 case does converge 
eventually, in practice we find the convergence is so slow that another Monte Carlo method 
is preferrable. Furthermore, the highest weights have the fewest events, and thus also suffer 
from small-number statistics, potentially adding to the "undue influence" of outliers. This 
can be seen in the scatter at the high- weight end of each data set. 

Instead of picking a scattering angle at random and weighting it by equation ( 147|) . let 
us use the differential cross section to get the scattering pdf: 

P(e,(f)) = —[{l-6){cos^e + l) + 26cosHcos^(f) + 26sm^(j)]. (Cll) 
IGtt 

Integrating over 0, we again find the standard Thomson cross section (this holds even for 

P{cos9) = l{cos^9 + l). (C12) 
8 

Writing /i = cos^ for convenience, the cumulative distribution function is given by 

cdf (/i) = r P{iJ,')dij,' = i(/i=^ + 3/i + 4) (C13) 
J-i 8 

To pick an appropriate value for fi, generate a random number uniformly distributed A G 
[0, 1), and invert equation flC13p . in effect solving for the root of the cubic: 

+ 3/i + 4 - 8A = 0. (C14) 
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Fig. 14. — Relative contribution to the total observed spectrum by photons of a given weight, 
for optical depths of r = (1, 2, 3, 5, 10). Any calculation with r > 2 will not formally converge 
using the method of weights. 
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Fig. 15. — Relative contribution to the total observed spectrum by photons of a given weight, 
for optical depths of r = (2, 5, 10), sorted by color (red, blue, black). The solid lines are the 
analytic results, and the crosses are "data" from Monte Carlo calculations of 10^ photons 
each. 
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Because cdf (/i) is monotonically increasing, this equation is guaranteed to have a single real 
root in the interval — 1 < < 1. 

Once 11 is selected, we choose by the same method, now using the pdf 

P(</); ii) = ^ , ^ , ^. [(1 - 5) (/i^ + 1) + 25/^2 cos^ + 25 sin^ 0] , (C15) 



cdf (0; ^) = ^ + —[ ^^^—^ ] sin(2</.). (C16) 



which gives 

27T ' A-n V/i^ + I, 

Again, to pick an appropriate (f) given a uniform random A, one must invert equation flC16l) 
to get = cdf~^(A). Unfortunately, this is equivalent to solving Kepler's equation, which 
has no closed-form solution, and must be done numerically. Fortunately, this is equivalent 
to solving Kepler's equation, one of the best-studied numerical problems i n astr ophysics! 



In practice, we use the iterative approach outlined in [Murray fc DermottI (jl999l ). While 
slightly more time consuming than the method of weights, the exact cross section method 
has the distinct advantage of converging for an arbitrary number of scatterings, and thus is 
the method we prefer for Pandurata. 
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